# Set up environment ----
library(data.table)
library(anesrake)
library(Hmisc)
library(ggplot2)
library(cowplot)
library(AER)
library(survey)
library(mediation)
library(parallel)
options(mc.cores = 6)
source("functions.R")

# Data ----
analysis_data <- readRDS("analysis-data-afpol.RDS")
analysis_data_uu <- readRDS("analysis-data-afpol-uninterested-unavailable.RDS")
ces2022 <- readRDS("ces2022pid7only.RDS")
data_for_flow_chart <- readRDS("data_for_flow_chart.RDS")

# Weight construction ----
targets <- list(
  region = c(
    Midwest = 0.2161336,
    Northeast = 0.1769097,
    South = 0.3798531,
    West = 0.2271037
  ),
  age_cat = c(
    `18-29` = 0.2088985,
    `30-44` = 0.2418020,
    `45-64` = 0.3322171,
    `65+` = 0.2170824
  ),
  education = c(
    `HS or Less` = 0.3705373,
    `Some College` = 0.3152755,
    `College Grad` = 0.1985126,
    `Post-Grad` = 0.1156746
  ),
  gender = c(
    Woman = 0.518939808691599,
    Man = 0.473758199179656,
    `other gender` = 0.00730199212874556
  ),
  race = c(
    `All Other` = 0.090643382,
    Asian = 0.046452363,
    Black = 0.126114272,
    `Native American` = 0.007927038,
    White = 0.728862944
  ),
  hispanic = c(
    Yes = 0.1395287,
    No = 0.8604713
  )
)
rake_fit <- rake(targets, analysis_data, analysis_data[, id])
analysis_data[, weight := rake_fit$weightvec]
des <- svydesign(~ id, data = analysis_data, weights = ~ weight)
rake_fit_uu <- rake(targets, analysis_data_uu, analysis_data_uu[, id])
analysis_data_uu[, weight := rake_fit_uu$weightvec]
des_uu <- svydesign(~ id, data = analysis_data_uu, weights = ~ weight)
design_list <- list(
  subset(des, !is.na(affective_polarization_voters)),
  subset(des, !is.na(in_party_affect_voters)),
  subset(des, !is.na(out_party_affect_voters)),
  subset(des, !is.na(affective_polarization_elites)),
  subset(des, !is.na(in_party_affect_elites)),
  subset(des, !is.na(out_party_affect_elites)),
  subset(des_uu, !is.na(affective_polarization_voters)),
  subset(des_uu, !is.na(in_party_affect_voters)),
  subset(des_uu, !is.na(out_party_affect_voters)),
  subset(des_uu, !is.na(affective_polarization_elites)),
  subset(des_uu, !is.na(in_party_affect_elites)),
  subset(des_uu, !is.na(out_party_affect_elites))
)

# Balance tests ----
all_X <- cbind(
  as.data.table(model.matrix(~ 0 + region, analysis_data)),
  as.data.table(model.matrix(~ 0 + age_cat, analysis_data)),
  as.data.table(model.matrix(~ 0 + hispanic, analysis_data)),
  as.data.table(model.matrix(~ 0 + race, analysis_data)),
  as.data.table(model.matrix(~ 0 + gender, analysis_data)),
  as.data.table(model.matrix(~ 0 + education, analysis_data)),
  as.data.table(model.matrix(~ 0 + work_status, analysis_data)),
  as.data.table(model.matrix(~ 0 + snap_status, analysis_data)),
  as.data.table(model.matrix(~ 0 + home_ownership, analysis_data)),
  as.data.table(model.matrix(~ 0 + political_interest, analysis_data)),
  as.data.table(model.matrix(~ 0 + partyid, analysis_data)),
  as.data.table(model.matrix(~ 0 + ideology, analysis_data)),
  pre_congress_approval = analysis_data$pre_congress_approval,
  pre_congress_trust = analysis_data$pre_congress_trust
)
missingness_vars <- which(colnames(all_X) %like% "missing|Independent")
all_X <- all_X[, -missingness_vars, with = FALSE]
weighted_equivalence_tests <- run_equiv(
  X = all_X,
  Tr = analysis_data$invited,
  w = analysis_data$weight,
  epsilon.method = "strict",
  fdr_correct = TRUE
)
unweighted_equivalence_tests <- run_equiv(
  X = all_X,
  Tr = analysis_data$invited,
  w = rep(1, nrow(analysis_data)),
  epsilon.method = "strict",
  fdr_correct = TRUE
)

# Formulas ----
lm_formulas <- list(
  affective_polarization_voters ~
    invited + org,
  in_party_affect_voters ~
    invited + org,
  out_party_affect_voters ~
    invited + org,
  affective_polarization_elites ~
    invited + org,
  in_party_affect_elites ~
    invited + org,
  out_party_affect_elites ~
    invited + org,
  affective_polarization_voters ~
    invited + org + interested + available,
  in_party_affect_voters ~
    invited + org + interested + available,
  out_party_affect_voters ~
    invited + org + interested + available,
  affective_polarization_elites ~
    invited + org + interested + available,
  in_party_affect_elites ~
    invited + org + interested + available,
  out_party_affect_elites ~
    invited + org + interested + available
)
iv_formulas <- list(
  affective_polarization_voters ~
    attended + org |
    invited + org,
  in_party_affect_voters ~
    attended + org |
    invited + org,
  out_party_affect_voters ~
    attended + org |
    invited + org,
  affective_polarization_elites ~
    attended + org |
    invited + org,
  in_party_affect_elites ~
    attended + org |
    invited + org,
  out_party_affect_elites ~
    attended + org |
    invited + org,
  affective_polarization_voters ~
    attended + org + interested + available |
    invited  + org + interested + available,
  in_party_affect_voters ~
    attended + org + interested + available |
    invited  + org + interested + available,
  out_party_affect_voters ~
    attended + org + interested + available |
    invited  + org + interested + available,
  affective_polarization_elites ~
    attended + org + interested + available |
    invited  + org + interested + available,
  in_party_affect_elites ~
    attended + org + interested + available |
    invited  + org + interested + available,
  out_party_affect_elites ~
    attended + org + interested + available |
    invited  + org + interested + available
)
ua_formulas <- list(
  affective_polarization_voters ~
    attended |
    invited,
  in_party_affect_voters ~
    attended |
    invited,
  out_party_affect_voters ~
    attended |
    invited,
  affective_polarization_elites ~
    attended |
    invited,
  in_party_affect_elites ~
    attended |
    invited,
  out_party_affect_elites ~
    attended |
    invited
)
sr_formulas <- list(
  affective_polarization_voters ~
    self_report_attended + org |
    invited + org,
  in_party_affect_voters ~
    self_report_attended + org |
    invited + org,
  out_party_affect_voters ~
    self_report_attended + org |
    invited + org,
  affective_polarization_elites ~
    self_report_attended + org |
    invited + org,
  in_party_affect_elites ~
    self_report_attended + org |
    invited + org,
  out_party_affect_elites ~
    self_report_attended + org |
    invited + org
)
dem_formulas <- list(
  affective_polarization_voters ~
    attended * dem + org |
    invited * dem + org,
  in_party_affect_voters ~
    attended * dem + org |
    invited * dem + org,
  out_party_affect_voters ~
    attended * dem + org |
    invited * dem + org,
  affective_polarization_elites ~
    attended * dem + org |
    invited * dem + org,
  in_party_affect_elites ~
    attended * dem + org |
    invited * dem + org,
  out_party_affect_elites ~
    attended * dem + org |
    invited * dem + org
)
strong_formulas <- list(
  affective_polarization_voters ~
    attended * strong + org |
    invited * strong + org,
  in_party_affect_voters ~
    attended * strong + org |
    invited * strong + org,
  out_party_affect_voters ~
    attended * strong + org |
    invited * strong + org,
  affective_polarization_elites ~
    attended * strong + org |
    invited * strong + org,
  in_party_affect_elites ~
    attended * strong + org |
    invited * strong + org,
  out_party_affect_elites ~
    attended * strong + org |
    invited * strong + org
)
mediation_base_formulas <- list(
  scaled_days_to_complete ~
    invited * org,
  scaled_days_to_complete ~
    invited * org,
  scaled_days_to_complete ~
    invited * org,
  scaled_days_to_complete ~
    invited * org,
  scaled_days_to_complete ~
    invited * org,
  scaled_days_to_complete ~
    invited * org
)
mediation_final_formulas <- list(
  affective_polarization_voters ~
    invited * org + scaled_days_to_complete * org,
  in_party_affect_voters ~
    invited * org + scaled_days_to_complete * org,
  out_party_affect_voters ~
    invited * org + scaled_days_to_complete * org,
  affective_polarization_elites ~
    invited * org + scaled_days_to_complete * org,
  in_party_affect_elites ~
    invited * org + scaled_days_to_complete * org,
  out_party_affect_elites ~
    invited * org + scaled_days_to_complete * org
)

# Estimation ----
linear_models <- lapply(1:12, function(x)
    summary(svyglm(lm_formulas[[x]], design = design_list[[x]]))
)
ivreg_models <- c(
  lapply(1:6, function(x)
    summary(svyivreg(ua_formulas[[x]], design = design_list[[x]]))),
  lapply(1:6, function(x)
    summary(AER::ivreg(iv_formulas[[x]], data = analysis_data))),
  lapply(1:12, function(x)
    summary(svyivreg(iv_formulas[[x]], design = design_list[[x]])))
)
slfrpt_models <- lapply(1:6, function(x)
    summary(svyivreg(sr_formulas[[x]], design = design_list[[x]]))
)
dem_models <- lapply(1:6, function(x)
    summary(svyivreg(dem_formulas[[x]], design = design_list[[x]]))
)
strong_models <- lapply(1:6, function(x)
    summary(svyivreg(strong_formulas[[x]], design = design_list[[x]]))
)
lm_days_to_complete_model <- svyglm(
  days_to_complete ~ invited + org, design = des
)
mediation_base_models <- lapply(1:6, function(x)
  svyglm(mediation_base_formulas[[x]], design = design_list[[x]])
)
mediation_final_models <- lapply(1:6, function(x)
  svyglm(mediation_final_formulas[[x]], design = design_list[[x]])
)
set.seed(2138051988, "L'Ecuyer-CMRG")
mediation_fit_list <- parallel::mclapply(1:6, function(i)
  mediate(
    mediation_base_models[[i]],
    mediation_final_models[[i]],
    sims = 1e3,
    treat = "invited",
    mediator = "scaled_days_to_complete"
  )
)

# Collect results ----
outcomes <- c(
  "Thermometer Difference (Voters)",
  "Inparty Warmth (Voters)",
  "Outparty Warmth (Voters)",
  "Thermometer Difference (Elites)",
  "Inparty Warmth (Elites)",
  "Outparty Warmth (Elites)"
)
estimands <- c(
  "Intent to Treat",
  "Complier Average Causal Effect",
  "CACE (Self reported attendance)"
)
adjusted_levels <- c("Unadjusted", "Adjusted")
weighted_levels <- c("Unweighted", "Weighted")
sample_levels <- c(
  "Exclude Uninterested or Unavailable",
  "Include Uninterested or Unavailable"
)
results <- CJ(
  estimand = factor(estimands, levels = estimands),
  adjusted = factor(adjusted_levels, levels = adjusted_levels),
  weighted = factor(weighted_levels, levels = weighted_levels),
  sample = factor(sample_levels, levels = sample_levels),
  outcome = factor(outcomes, levels = rev(outcomes)),
  sorted = FALSE)[
    !(estimand %like% "^I" & (adjusted %like% "U" | weighted %like% "U")) &
    !(estimand %like% "^CA" & (adjusted %like% "U" | weighted %like% "U")) &
    !(adjusted %like% "U" & weighted %like% "U") &
    !(sample %like% "I" & adjusted %like% "U") &
    !(sample %like% "I" & weighted %like% "U") &
    !(sample %like% "I" & !estimand %like% "Co|I")
  ]
results[estimand %like% "I", `:=`(
  estimate = sapply(linear_models, function(x)
    x$coef["invited", "Estimate"]),
  std_error = sapply(linear_models, function(x)
    x$coef["invited", "Std. Error"]),
  pvalue = sapply(linear_models, function(x)
    x$coef["invited", "Pr(>|t|)"])
)]
results[estimand %like% "E" & estimand %like% "Co", `:=`(
  estimate = sapply(ivreg_models, function(x)
    x$coef["attended", "Estimate"]),
  std_error = sapply(ivreg_models, function(x)
    x$coef["attended", "Std. Error"]),
  pvalue = sapply(ivreg_models, function(x)
    x$coef["attended", "Pr(>|t|)"])
)]
results[estimand %like% "E" & estimand %like% "CA", `:=`(
  estimate = sapply(slfrpt_models, function(x)
    x$coef["self_report_attended", "Estimate"]),
  std_error = sapply(slfrpt_models, function(x)
    x$coef["self_report_attended", "Std. Error"]),
  pvalue = sapply(slfrpt_models, function(x)
    x$coef["self_report_attended", "Pr(>|t|)"])
)]
denom_formulas <- list(
  ~ affective_polarization_voters,
  ~ in_party_affect_voters,
  ~ out_party_affect_voters,
  ~ affective_polarization_elites,
  ~ in_party_affect_elites,
  ~ out_party_affect_elites
)
results[weighted %like% "W" & sample %like% "E",
  cohens_d_denom := rep(sapply(1:6, function(x) {
    svyvar(denom_formulas[[x]], design_list[[x]]) ^ .5
  }), 4)]
results[weighted %like% "W" & sample %like% "I",
  cohens_d_denom := rep(sapply(1:6, function(x) {
    svyvar(denom_formulas[[x]], design_list[[x]]) ^ .5
  }), 2)]
results[weighted %like% "U",
  cohens_d_denom := c(
    sd(analysis_data$affective_polarization_voters, na.rm = TRUE),
    sd(analysis_data$in_party_affect_voters, na.rm = TRUE),
    sd(analysis_data$out_party_affect_voters, na.rm = TRUE),
    sd(analysis_data$affective_polarization_elites, na.rm = TRUE),
    sd(analysis_data$in_party_affect_elites, na.rm = TRUE),
    sd(analysis_data$out_party_affect_elites, na.rm = TRUE)
  )]
results[, cohens_d := abs(estimate / cohens_d_denom)]
results[, outcome := factor(outcome, levels = outcomes)]
results[, fdr_p := p.adjust(pvalue, "fdr"),
  .(estimand, adjusted, weighted, sample)]
results[, pretty_fdr_p := sprintf("%.04f", round(fdr_p, 4))]
results[, pretty_cohens_d := sprintf("%.02f", round(cohens_d, 2))]
results[, pretty_label := sprintf("%.0f", estimate)]
results[, pretty_pvalue := sprintf("%.04f", round(pvalue, 4))]

# Generate figures ----
ggdata <- copy(results)
ggdata[, outcome_name := gsub("\\(Elites\\)|\\(Voters\\)", "", outcome)]
ggdata[, object_name := ifelse(outcome %like% "Elite", "Elites", "Voters")]
ggdata[, self_reported := factor(ifelse(estimand %like% "CACE",
  "Self-Reported Attendance", "Verified Attendance"))]
ggdata[estimand == "CACE (Self reported attendance)",
  estimand := "Complier Average Causal Effect"]
ggdata[, group := paste(adjusted, weighted, sample, self_reported, sep = "/")]
v <- 1
shape_values <- c(19, 15, 17, 23, 25)
names(shape_values) <- c(
  paste0("Adjusted/Weighted/Exclude Uninterested or Unavailable/",
    "Verified Attendance"),
  paste0("Adjusted/Weighted/Include Uninterested or Unavailable/",
    "Verified Attendance"),
  paste0("Unadjusted/Weighted/Exclude Uninterested or Unavailable/",
    "Verified Attendance"),
  paste0("Adjusted/Unweighted/Exclude Uninterested or Unavailable/",
    "Verified Attendance"),
  paste0("Adjusted/Weighted/Exclude Uninterested or Unavailable/",
    "Self-Reported Attendance")
)
ggplot(ggdata,
  aes(outcome_name, estimate, group = group, shape = group)) +
  geom_hline(yintercept = 0, linetype = 3, color = "lightgray") +
  geom_errorbar(aes(
    ymin = estimate + qnorm(.025) * std_error,
    ymax = estimate + qnorm(.975) * std_error),
    position = position_dodge(width = v), width = 0) +
  geom_point(position = position_dodge(width = v), size = 2, fill = "black") +
  coord_flip() +
  ylab("Estimate (95% Interval)") +
  xlab("") +
  facet_grid(object_name ~ estimand, scales = "free_x") +
  theme_bw() +
  ggtitle("") +
  scale_shape_manual(values = shape_values) +
  theme(
    strip.background = element_rect(fill = NA, color = NA),
    legend.position = "bottom",
    plot.title = element_text(hjust = .5),
    legend.title = element_blank(),
    legend.direction = "vertical")
ggsave("figure-afpol-all.png", width = 6.5, height = 5)

v <- .4
ggplot(ggdata[weighted %like% "W" & adjusted %like% "A" &
    self_reported %like% "V" & sample %like% "E"],
  aes(outcome_name, estimate, group = group)) +
  geom_hline(yintercept = 0, linetype = 3, color = "lightgray") +
  geom_errorbar(aes(
    ymin = estimate + qnorm(.025) * std_error,
    ymax = estimate + qnorm(.975) * std_error),
    position = position_dodge(width = v), width = 0) +
  geom_point(position = position_dodge(width = v), size = 2) +
  coord_flip() +
  ylab("Estimate (95% Interval)") +
  xlab("") +
  facet_grid(object_name ~ estimand, scales = "free_x") +
  theme_bw() +
  ggtitle("") +
  theme(
    strip.background = element_rect(fill = NA, color = NA),
    legend.position = "bottom",
    plot.title = element_text(hjust = .5),
    legend.title = element_blank()) +
  geom_label(aes(label = pretty_label))
ggsave("figure-afpol.png", width = 6.5, height = 5)

dn <- c("Midwest", "Northeast", "South", "West", "18-29", "30-44", "45-64",
  "65+", "Hispanic/Latino", "Not Hispanic/Latino", "All Other Races/Mixed",
  "Asian", "Black", "Native American", "White", "Woman", "Man", "Other gender",
  "High School Grad or Less", "Some College", "College Grad", "Post-Grad",
  "Works full-time", "Works part-time", "Does not work", "Receives SNAP",
  "Does not receive SNAP", "Home owned with loan", "Home owned without loan",
  "Home rented", "Home occupied without payment",
  "Follows politics most of the time", "Follows politics some of the time",
  "Follows politics only now and then", "Follows politics hardly at all",
  "Strong Democrat", "Not Strong Democrat", "Lean Democrat", "Lean Republican",
  "Not Strong Republican", "Strong Republican", "Conservative",
  "Extremely conservative", "Extremely liberal",
  "Haven't thought much about this", "Liberal", "Moderate; middle of the road",
  "None of these", "Slightly conservative", "Slightly liberal",
  "Congress approval (pre-event)", "Congress trust (pre-event)")
dn_ordering <- c(5:8, 19:22, 16:18, 9:10, 11:15, 1:4, 23:31, 32:41, 43, 42, 49,
  47, 50, 46, 44, 45, 48, 51:52)
generate_equivalence_test_plot(
  weighted_equivalence_tests[dn_ordering],
  display_names = dn[dn_ordering],
  panel_widths = c(1, .4, 2, .4, .4),
  title_text = "Equivalence Tests of Balance (Weighted)")
ggsave("figure-afpol-balance-weighted.png", width = 6.5, height = 7,
  scale = 1.67)

generate_equivalence_test_plot(
  unweighted_equivalence_tests[dn_ordering], display_names = dn[dn_ordering],
  panel_widths = c(1, .4, 2, .4, .4),
  title_text = "Equivalence Tests of Balance (Unweighted)")
ggsave("figure-afpol-balance-unweighted.png", width = 6.5, height = 7,
  scale = 1.67)

# Quantities reported in main text ----
# Respondents invited to the DTH...
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Elites)", pretty_label]
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Voters)", pretty_label]
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Elites)", pretty_cohens_d]
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Voters)", pretty_cohens_d]
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Elites)", pretty_fdr_p]
results[estimand %like% "I" & sample %like% "E" &
    outcome == "Thermometer Difference (Voters)", pretty_fdr_p]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" & outcome == "Thermometer Difference (Elites)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" &outcome == "Thermometer Difference (Voters)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]

# To better understand these results, ...
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" &outcome == "Outparty Warmth (Elites)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" &outcome == "Inparty Warmth (Elites)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" & outcome == "Outparty Warmth (Voters)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" & outcome == "Inparty Warmth (Voters)", .(
      pretty_label, pretty_cohens_d, pretty_fdr_p)]

# Results are substantively robust...
results[estimand %like% "CA" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" &outcome %like% "Thermometer Difference \\(E", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "CA" & adjusted %like% "A" & weighted %like% "W" &
    sample %like% "E" &outcome %like% "Thermometer Difference \\(V", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "U" &
    sample %like% "E" & outcome %like% "Thermometer Difference \\(E", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]
results[estimand %like% "Co" & adjusted %like% "A" & weighted %like% "U" &
    sample %like% "E" & outcome %like% "Thermometer Difference \\(V", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]
results[estimand %like% "Co" & adjusted %like% "U" & weighted %like% "W" &
    sample %like% "E" & outcome %like% "Thermometer Difference \\(E", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]
results[estimand %like% "Co" & adjusted %like% "U" & weighted %like% "W" &
    sample %like% "E" & outcome %like% "Thermometer Difference \\(V", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]
results[estimand %like% "Co" & sample %like% "I" &
    outcome %like% "Thermometer Difference \\(E", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]
results[estimand %like% "Co" & sample %like% "I" &
    outcome %like% "Thermometer Difference \\(V", .(
      outcome, pretty_label, pretty_cohens_d, pretty_fdr_p, pretty_pvalue)]

# The failure to reject the null for the unweighted analysis...
round(min(p.adjust(sapply(dem_models,
  function(x) x$coef["attended:dem", 4]), "fdr")), 4)
round(min(sapply(dem_models,
  function(x) x$coef["attended:dem", 4])), 4)
round(min(p.adjust(sapply(strong_models,
  function(x) x$coef["attended:strong", 4]), "fdr")), 4)
round(min(sapply(strong_models,
  function(x) x$coef["attended:strong", 4])), 4)

# Furthermore, we observe no evidence consistent...
round(summary(lm_days_to_complete_model)$coef["invited", -3], 2)
round(range(sapply(mediation_fit_list, function(x) summary(x)$d1)), 2)
min(p.adjust(sapply(mediation_fit_list, function(x) summary(x)$d1.p), "fdr"))
min(sapply(mediation_fit_list, function(x) summary(x)$d1.p))

# Supporting Materials ----
#   Sample description ----
n <- analysis_data[, .N]
cat(
  "We recruited samples from three survey houses: CHRR at Ohio State Universit",
  "y, SurveyUSA, and YouGov. For the sample in the main analyses, the demograp",
  "hic and geographic breakdown is: [Age] ",
  analysis_data[, round(100 * .N / n, 0), age_cat][order(age_cat)][,
    paste0(paste0(age_cat, " (", V1, "%)"), collapse = ", ")], "; [Education] ",
  analysis_data[, round(100 * .N / n, 0), education][order(education)][,
    paste0(paste0(education, " (", V1, "%)"), collapse = ", ")], "; [Gender] ",
  analysis_data[, round(100 * .N / n, 0), gender][order(gender)][,
    paste0(paste0(gender, " (", V1, "%)"), collapse = ", ")],
  "; [Hispanic/Latino status] ",
  analysis_data[, round(100 * .N / n, 0), hispanic][order(hispanic)][,
    paste0(paste0(hispanic, " (", V1, "%)"), collapse = ", ")], "; [Race] ",
  analysis_data[, round(100 * .N / n, 0), race][order(race)][,
    paste0(paste0(race, " (", V1, "%)"), collapse = ", ")], "; [Region]",
  analysis_data[, round(100 * .N / n, 0), region][order(region)][,
    paste0(paste0(region, " (", V1, "%)"), collapse = ", ")], ".\n\n",

  "The breakdown by political variables is [Party ID] ",
  analysis_data[, round(100 * .N / n, 0), partyid][order(partyid)][,
    paste0(paste0(partyid, " (", V1, "%)"), collapse = ", ")],
  "; [Political Ideology] ",
  analysis_data[, round(100 * .N / n, 0), ideology][
    c(6, 1, 4, 2, 7, 3, 5, 8:9)][,
    paste0(paste0(ideology, " (", V1, "%)"), collapse = ", ")],
  "; [Political Interest] ",
  analysis_data[, round(100 * .N / n, 0), political_interest][c(1:3, 5)][,
      paste0(paste0(political_interest, " (", V1, "%)"), collapse = ", ")],
  ".",
  sep = "")

#   Weight construction ----
wa <- anesrake::weightassess(targets, analysis_data, rake_fit$weightvec)
pretty_wa <- lapply(wa, function(x) {
  paste0(paste0(
    rownames(head(x, -1)), " (",
    gsub("-", "\u2212", round(100 * head(x, -1)[, 1], 0)), "%; ",
    gsub("-", "\u2212", round(100 * head(x, -1)[, 8], 0)), "%)"),
    collapse = ", ")
})
ces_denom <- ces2022[pid7 %like% "D|R", sum(commonweight)]
party_line <- cbind(
  analysis_data[,
    round(100 * .N / n, 0), partyid][order(partyid)],
  V2 = analysis_data[,
    round(100 * sum(weight) / n, 0), partyid][order(partyid)]$V1,
  V3 = ces2022[pid7 %like% "D|R",
    round(100 * sum(commonweight) / ces_denom, 0), pid7][order(pid7)]$V1)[,
      paste0(paste0(partyid,
        c(" (unweighted = ", rep(" (", 5)), V1, "%, ",
        c("weighted = ", rep("", 5)), V2, "%, ",
        c("target = ", rep("", 5)), V3, "%",
        ")"), collapse = ", ")]
cat(
  "We targeted age category, educational attainment, gender, Hispanic/Latino s",
  "tatus, race, and census region. The targets and discrepancies (i.e., target",
  " – actual) for these variables were [Age] ",
  pretty_wa$age_cat, "; [Education] ",
  pretty_wa$education, "; [Gender] ",
  pretty_wa$gender, "; [Hispanic/Latino Status] ",
  pretty_wa$hispanic, "; [Race] ",
  pretty_wa$race, "; [Region] ",
  pretty_wa$region, ". ",
  "Weights range from ",
  round(min(rake_fit$weightvec), 2), " to ",
  round(max(rake_fit$weightvec), 2), ", with a median of ",
  round(median(rake_fit$weightvec), 2), " and a mean of ",
  round(mean(rake_fit$weightvec), 2), ". ",
  "Weighting by these demographic and geographic covariates also improved repr",
  "esentativeness on Party ID for strong identifiers, with unweighted margins,",
  "weighted margins, and targets from the 2022 Cooperative Election Study resp",
  "ectively as follows: ",
  party_line,
  ".",
  sep = "")

#   Balance ----
weighted_equivalence_tests[, .N]
weighted_equivalence_tests[fdr_p < .05, .N]
unweighted_equivalence_tests[fdr_p < .05, .N]
weighted_equivalence_tests[, .(max(fdr_p, na.rm = TRUE), sum(!is.na(fdr_p)))]
unweighted_equivalence_tests[, .(max(fdr_p, na.rm = TRUE), sum(!is.na(fdr_p)))]

#   Self-reported attendance analysis ----
print_results(results[estimand %like% "CA"][c(4:6, 1:3)])

#   Unweighted analysis ----
print_results(results[weighted %like% "U"][c(4:6, 1:3)])

#   Unadjusted analysis ----
print_results(results[adjusted %like% "U"][c(4:6, 1:3)])

#   Including uninterested/unavailable ----
print_results(results[estimand %like% "I" & sample %like% "I"][c(4:6, 1:3)])
print_results(results[estimand %like% "Co" & sample %like% "I"][c(4:6, 1:3)])

#   Heterogeneity ----
lapply(c(4:6, 1:3), print_het_results, models = dem_models,
  var = "attended:dem")
round(min(p.adjust(sapply(dem_models,
  function(x) x$coef["attended:dem", 4]), "fdr")), 4)
lapply(c(4:6, 1:3), print_het_results, models = strong_models,
  var = "attended:strong")
round(min(p.adjust(sapply(strong_models,
  function(x) x$coef["attended:strong", 4]), "fdr")), 4)

#   Decay mediation ----
lapply(c(4:6, 1:3), print_med_results, models = mediation_fit_list)
round(min(p.adjust(sapply(mediation_fit_list,
  function(x) x$d1.p), "fdr")), 4)

#   Generate table ----
results[, .(
  outcome,
  estimand,
  adj = factor(ifelse(adjusted == "Adjusted", "Yes", "No"),
    levels = c("Yes", "No")),
  wtd = factor(ifelse(weighted == "Weighted", "Yes", "No"),
    levels = c("Yes", "No")),
  int_av = factor(ifelse(sample %like% "I", "Yes", "No"),
    levels = c("Yes", "No")),
  est = pretty_label,
  se = round(std_error, 0),
  f = pretty_cohens_d,
  pfdr = pretty_fdr_p,
  uncp = sprintf("%.04f", pvalue))][order(outcome, estimand, adj, wtd)]

#   Flow chart (data for Fig 1, figure created in PowerPoint) ----
data_for_flow_chart[, .N, .(interested, available, continue_with_study, invited,
  in_post, attended, asked_ft)][order(interested, available,
    continue_with_study, invited, in_post, attended, asked_ft)]
